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Abstract. We study the mean- field dynamics and the reduced-dimension character 
of two-mode Bose-Einstein condensates (BECs) in highly anisotropic traps. By means 
of perturbative techniques, we show that the tightly confined (transverse) degrees of 
freedom can be decoupled from the dynamical equations at the expense of introducing 
additional effective three-body, attractive, intra- and inter-mode interactions into the 
dynamics of the loosely confined (longitudinal) degrees of freedom. These effective 
interactions are mediated by changes in the transverse wave function. The perturbation 
theory is valid as long as the nonlinear scattering energy is small compared to 
the transverse energy scales. This approach leads to reduced-dimension mean-field 
equations that optimally describe the evolution of a two-mode condensate in general 
quasi-lD and quasi-2D geometries. We use this model to investigate the relative 
phase and density dynamics of a two- mode, cigar-shaped ^''Rb BEC. We study the 
relative- phase dynamics in the context of a nonlinear Ramsey interferometry scheme, 
which has recently been proposed as a novel platform for high-precision interferometry. 
Numerical integration of the coupled, time-dependent, three-dimensional, two-mode 
Gross-Pitaevskii equations for various atom numbers shows that this model gives 
a considerably more refined analytical account of the mean-field evolution than an 
idealized quasi- ID description. 
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1. Introduction 

Quantum protocols for nonlinear interferometry using two-mode Bose-Einstein 
condensates (BECs) have recently been proposed as a novel platform for weak-signal 
detection [1, 2, 3]. More specifically a two-mode ^''Rb BEC consisting of atoms 
can be used to implement a nonlinear Ramsey interferometer that potentially operates 
near the limits established by quantum mechanics [4, 5, 6]. This protocol is able to 
achieve detection sensitivities that scale better than the optimal 1/A^ limit of linear 
interferometry [ ], without relying on complicated state preparation or measurement 
procedures nor on entanglement generation to enhance the measurement sensitivity [ ] . 
For these reasons, this scheme can be particularly attractive from an experimental 
perspective. 

As discussed in previous work [4, 5, 6], there are several conditions for observing 
nonlinearity-enhanced scalings in this system. For instance, the expansion of the 
condensate with increasing atom number can get in the way of achieving the desired 
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scaling, as it essentially dilutes the nonlinear interactions in the condensate. Although 
this effect can be completely suppressed by hard- walled homogeneous potentials, 
confining the condensate in highly anisotropic traps is a more realistic solution. In this 
way, the expansion of the atomic cloud in the tightly confined directions is effectively 
suppressed, with highly elongated geometries providing the architecture with best 
scalings. 

In view of currently available techniques and typical experimental parameters, we 
numerically simulated such a nonlinear BEC interferometer in cigar-shaped (quasi-lD) 
potentials [ ]. Interestingly, our simulations revealed that the interferometer's signal can 
be sensitive to the three-dimensional nature of the condensate. In spite of the highly 
elongated confinement, we found significant deviations from a quasi- ID model as the 
strength of the nonlinear scattering interaction increases. In this reduced-dimension 
approximation, it is assumed that the tightly confined (transverse) dimensions can be 
effectively neglected on the grounds that the characteristic transverse energy scale far 
exceeds the scattering interaction energy of the atomic cloud. In this situation, however, 
there are still position-dependent phase shifts that need to be modeled precisely for an 
accurate analytical description of the interferometry process. This brings into question 
the accuracy of the reduced-dimension approximation, both spatially and temporally. 

To better model the ground-state properties of a BEC in highly anisotropic 
potentials, we derived corrections to such approximation by means of perturbative 
techniques [ ]. Using a perspective borrowed from quantum information theory, we 
developed a perturbative Schmidt decomposition of the condensate wave function 
between the transverse direction(s) and the loosely confined (longitudinal) direction(s). 
This formalism provides corrections to the lowest-order transverse and longitudinal wave 
functions of the reduced-dimension approximation; the main effect is a reshaping of the 
BEC in the tightly confined direction(s) as the strength of the nonlinear scattering 
interaction increases. In addition, because the perturbation formalism is tied to the 
Schmidt decomposition, it automatically encodes information about the entanglement 
between the spatial directions in higher-order Schmidt terms; the leading Schmidt term 
provides the optimal product approximation to the exact three-dimensional ground-state 
mean-field solution. 

In this article, we develop further the approach of Ref. [ ] by extending our 
perturbation theory to the mean-field dynamics of two-mode BECs in highly anisotropic 
traps. The perturbation theory is valid as long as the nonlinear scattering energy is 
small compared to the transverse energy scales. This approach leads to equations that 
effectively describe the evolution of a two-mode condensate in general quasi- ID and 
quasi-2D geometries. In addition, we also derive effective time-evolution equations for a 
highly anisotropic, single-mode BEC. These equations show how the corrections to the 
reduced-dimension approximation propagate in time and affect the overall dynamics of 
the condensate, thus modeling not only interference effects, but also density dynamics. 
We apply this 3D-corrected model to analyze the nonlinear BEC interferometry protocol 
simulated in [ ] for cigar-shaped potentials. This analysis leads to a considerably more 
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refined analytical account of tiie interference signal than the one given by a quasi- 
ID model. In addition, we study the spatial segregation of the two modes, which 
occurs on a longer timescale than the phase accumulation. We compare this model 
against exact three-dimensional numerical results for the evolution of the two-mode 
BEC given by the three-dimensional, two-mode, coupled Gross-Pit aevskii equations. 
Finally, we investigate the reduced-dimension structure of the highly anisotropic, two- 
mode condensate by explicitly deriving the instantaneous Schmidt decomposition of 
the time-dependent, two-mode condensate state. This analysis allows us to study the 
entanglement between the transverse and the longitudinal-internal degrees of freedom. 
It confirms that the tightly confined dimensions can indeed be decoupled from the 
evolution equation and naturally shows the optimal way to do it. 

This article is organized as follows. In Sec. 2, we briefly recount the reduced- 
dimension approximation to the three-dimensional, two-mode, coupled Gross-Pitaevskii 
equations, which completely neglects the effect of the nonlinear scattering interaction 
on the transverse degrees of freedom. By developing a perturbative relative-state 
decomposition of the time-dependent condensate mean fields, we derive corrections to 
this approximation in Sec. 3; these corrections act as effective three-body, attractive, 
intra- and inter-mode interactions. Because of the high anisotropy of the condensate, 
the evolution of the longitudinal modes in the perturbative regime takes place on a 
much longer time scale than the transverse modes. We address this separation of times 
scales in Sec. 4 and derive an adiabatic approximation to the longitudinal evolution 
equations of Sec. 3. In Sec. 5, we use this model to investigate the relative-phase and 
density dynamics of a two-mode, cigar-shaped ^''Rb BEC. We study the relative-phase 
dynamics in the context of a nonlinear Ramsey interferometry scheme. We compare 
our perturbative model against numerical integration of the full three-dimensional two- 
mode, coupled Gross-Pitaevskii equations for various atom numbers. Lastly, in Sec. 6, we 
investigate the reduced-dimension character of the two-mode condensate by explicitly 
deriving the instantaneous Schmidt decomposition of the time-dependent, two-mode 
condensate singe-particle state. Conclusions are given in Sec. 7. 

2. Two- mode dynamics in the reduced-dimension approximation 

In the mean-field approximation, one describes the evolution of a two-mode condensate 
at zero temperature by unity-normalized wave functions '0^(r,t) and '0^(r,t) for the 
two modes, which are determined by the two-mode, coupled, time-dependent, three- 
dimensional Gross-Pitaevskii (CP) equations 

pjj.a . fc2 . 

= ( - + F + ga^{N^ - l)\r? + 9a0{Nf> - 1)1^^1^)^". (1) 

Here A^i(2) is the number of atoms in mode 1 (2), 1/ is the external trapping potential, 
assumed to be the same for both modes, and a — 1^2 and /3 = 2, 1 label the two 
modes of the condensate, which can represent, for instance, two hyperfine states, |1) 
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and 1 2). We further assume that coUisions between the atoms are elastic, so that the 
only aUowed scattering processes are 1 1)1 1) |1)|1), |2)|2) |2)|2), and |1)|2) |1)|2), 
with scattering strengths gn^ 5^22 , and gi2^ where gai3 = Anh^a^fi/m is determined by 
the 5- wave scattering length and the atomic mass m. Note that by setting the 
inter-mode coupling constant gi2 to zero, one recovers the case of a single-mode BEC, 
which we address in more detail in Sec. 4.1. 

In the case of highly anisotropic potentials, the condensate is loosely trapped by a 
potential Vl(^) in d dimensions, referred to as longitudinal (L) dimensions, as opposed to 
the remaining D = 3 — d transverse degrees of freedom (T) , which are tightly confined 
in a potential Vt{p)- If the scattering interaction is sufficiently small compared to 
the transverse energy scale, one can neglect the effect of the nonlinear interaction on 
the atomic transverse degrees of freedom and hence approximate the condensate wave 
functions by the product ansatz 



where £"0 and Co(p) are, respectively, the ground-state energy and wave function of 
the bare transverse potential, Vt{p). The longitudinal wave functions (j)^{r^t) are the 
solutions of the rf-dimensional, longitudinal GP equations. 



= -T^^l + VL+gaaMNa " lW\^ + Qa^N p - 1)|0^|2 U", (3) 



which are found by plugging Eq. (2) into the two-mode GP equations (1) and projecting 
the result onto the subspace spanned by (^o- Note that in this reduced-dimension picture, 
the coupling constants are renormalized by the average inverse transverse cross section 
of the condensate, given by 



Within this approximation, the transverse and longitudinal degrees of freedom 
are decoupled. This approximation is only meaningful if the number of atoms in the 
condensate is small compared to an (upper) critical atom number A/r, defined as the 
number at which the nonlinear scattering energies become comparable to the transverse 
kinetic energy. As approaches A/r, one can no longer neglect the effects of the 
scattering interaction on the condensate transverse degrees of freedom; as a result, the 
product ansatz (2) is no longer a good approximation to the 3D wave function. Such 3D- 
induced effects are responsible not only for modifying the transverse and longitudinal 
wave functions, but also for creating correlations between the spatial directions. All 
these effects can be readily calculated in the perturbative regime where N is small 
compared to A/r- 

3. Two-mode dynamics via the relative-state decomposition 



t) = e-^'^^''%{p)r{r^t) 



(2) 






(4) 



In previous work [ ], we found perturbative corrections to the reduced-dimension 
approximation for the ground-state wave function of a single-mode condensate. These 
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corrections arise from the derivation of the Schmidt decomposition of the condensate 
wave function between the transverse and longitudinal degrees of freedom. To first-order 
this approach is equivalent to a first-order perturbative relative-state decomposition of 
the condensate wave function. The crucial difference between these two approaches is 
that the Schmidt decomposition assumes no prior knowledge of the basis elements used 
in the perturbation expansion, whereas in the relative-state method, the expansion is 
carried out relative to a fixed basis for the transverse degrees of freedom. For this reason, 
the expansion of the time-independent condensate wave function can be implemented 
in a simpler way if derived via the relative-state method rather than via the Schmidt 
decomposition. This motivates us to investigate the dynamics of the condensate from 
a similar perspective and to look for the perturbative relative-state decomposition of 
the time-dependent condensate mean field. We return to the Schmidt decomposition in 
Sec. 6 and show that a Schmidt decomposition of the three-dimensional, time-dependent 
mean-field solution can easily be retrieved from the relative-state decomposition. 

As before, we model the two-mode dynamics of the condensate by means of the 
mean-field approximation, according to which the mean field of each BEC mode satisfies 
the time-dependent, coupled, two-mode Gross-Pitaevskii equation (1), which we write 
here as 

ihr = {Ht + cHl + eg^^\r? + e^.^l^^l^)^^ (5) 

For simplicity of notation, we now use ijj^ — dijj^/dt^ g^a = {^a — ^)9aa^ and 
ga(3 — {N(3 — ^)ga(3' Here Ht{l) — — (^^/2m)V^(^) +Vt(l) is the transverse (longitudinal) 
single-particle Hamiltonian, and e is a formal perturbation parameter that is set equal 
to 1 at the end of the calculation. 

In the perturbative regime, we write the relative-state decomposition of the time- 
dependent condensate wave functions ip*^ as 

oo 

rip, r, t) = eo(p)^o t)+^J2 UpX{r, t\ (6) 

n=l 

where {^n} is the eigenbasis of the transverse Hamiltonian. The longitudinal wave 
functions are defined by the projection of the time-dependent mean- field solution ijj^ 
onto the transverse eigenfunctions 

Before we carry out the perturbative expansion of the time-dependent, two-mode 
GP equations, it is convenient to redefine the longitudinal wave functions relative to an 
interaction picture in which fast, trivial oscillations are removed. Because of the high 
anisotropy of the trapping potential, we expect the transverse bare trap energy to be 
the fastest time scale in the perturbative regime, and therefore we define 

(p«(t)^e*W^«(t). (7) 
Now we expand the two-mode GP equations (5) to second order in powers of e and 
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CO 

oo 

+ e'^a/5<P0 E(^olen) (<P?Cpr + (p?*(p^) . (8) 

n=l 

Hereafter, for brevity, we represent spatial integrals in terms of bra-ket inner products. 

Similarly, the projection of Eq. (5) onto Cn, ^ ^ 1, yields to lowest order in e the 
time-evolution equations for the longitudinal functions (p^, which read 

ihi^l = {Er, - Eo)<P^ + i^nieo) {gaaM' + ~9aM\') <Po + 0{e), U > 1. (9) 

One can integrate Eq. (9), obtaining the formal solution 

(10) 

By plugging Eq. (10) back into Eq. (8), we eliminate (p"(t) from the equation and 

get 

ihdp^ = e(HL + gaaTirWo? + Qa^riTWo?) <Po 
+ e'r^^f2^„a|cPor + ^a/3kor) 

ds {gaaW'ois)? + ~gaM{s)?) Grit - s)(p?(s) 

+ eWT~9o.p^>^o ds [m\<PoisW + ~9^a\^>Us)\') (cp?G^(t - s)cp?*(s) 

+ (p^*Gr(t-s)(p^(5)) 
+ e%"(t). (11) 

Here 
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is a temporal response function that comes from changes in the transverse wave function 
and that osciUates at the bare transverse eigenfrequencies, and 

oo 
n=l 

oo 
n=l 

oo 

+ ~9apV^ (cp^(pf (0)e^(^»-^°)*/'^ + (p^ (p^(0)e-^(^"-^°)*/^) (13) 

n=l 

is the term associated with the initial condensate wave function. 

Equations (10) and (11) are the final result of the relative-state method. Not 
surprisingly, cpg satisfies to lowest order the idealized reduced-dimension description of 
the time-dependent GP equations discussed in Sec. 2. The higher-order terms, on the 
other hand, introduce corrections to the reduced-dimension approximation that act as 
effective three-body, attractive, intra- and inter-mode interactions. These interactions 
are mediated by the changes in the transverse wave function and take place on the 
transverse time scales contained in the temporal response function Gt- 
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Figure 1. Temporal response function Gt, in units of ^ for a pancake-like trap 
{D = 1) and a transverse harmonic potential. Time is measured in units of uj^^ . 
The dashed, black hue represents the real part of Gt, the dotted, black hue is the 
imaginary part, and the solid, purple line is the absolute value of Gt- Gt is a periodic 
function with fundamental oscillation frequency equal to 2ujt- The anharmonicity 
of Gt essentially comes from the lowest two nonvanishing modes in the spectral 
decomposition (12), as shown in Fig. 3. 

Interestingly, the temporal response function (12) has a closed-form solution for the 
case of a transverse harmonic potential. In this case, the transverse ground-state wave 
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Figure 2. Temporal response function Gt, in units of ^ for a cigar-shaped trap 
{D = 2) and a transverse harmonic potential. Time is measured in units of cj^^. The 
dashed, blue line represents the real part of Gt, the dotted, blue line corresponds to the 
imaginary part, and the solid, orange line is the absolute value of Gt- As in the D = 1 
case, the temporal response function Gt for a cigar-shaped BEC is a periodic function 
with fundamental oscillation frequency equal to 2ujt- The anharmonicity of Gt is 
mainly due to the lowest two nonvanishing modes in the spectral decomposition (12) 
(see Fig. 3). 



function is the Gaussian 

where po = \/h/mujT^ with ujt being the transverse trap frequency. It is easy to see 
that 

Moreover, for a pancake-like potential {D = 1), we have 

{-IT" r/^^ + l^ 

-r]T^ , even, 



{in\il)={ ^\ " V 2 y ' (16) 

0, n odd. 

For a cigar (L) = 2), if we use polar coordinates for the transverse eigenf unctions, they 
take the form ^nrm{p^ 9^), with and m being radial and azimuthal quantum numbers 
and with the eigenenergies given by En^^ = hujT{2nr + \m\ + 1). Then we find that 

iUmieoo) = 2-"'-r]T5„zo- (17) 
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By using Eqs. (16) and (17), we can write the response function as 

2 



Grit) = { 



1 

1 



h 1 - Ae^i^Tt ' 



D = \ (pancake), 
D = 2 (cigar). 



(18) 




n 



Figure 3. Value of the coefficients (^ol^n)^ (in units of r]^) for a harmonic transverse 
trap for the first eight terms in the series (12). For every n, we show the value of the 
coefficient for a cigar-shaped potential (black bar on the left) as well as for a pancake- 
like trap (blue bar on the right). For a quasi- ID potential {D = 2), n represents 
the radial quantum number of the harmonic potential. For both cases, the spectral 
representation of Gt reveals a single dominant mode, namely n = 1 for a cigar and 
n = 2 for a pancake-like potential, that oscillates with frequency equal to 2cjt- The 
remaining, smaller terms introduce the small anharmonicity seen in Figs. 1 and 2. 

Figures 1 and 2 show the temporal response function (18) for 2D and ID trap 
geometries. For both trap geometries, the fundamental oscillation frequency of Gt 
is equal to 2(jJT' The anharmonicity of Gt is better understood from the spectral 
representation (12). Figure 3 shows the value of the coefficients (ColCn)^ for the first 
few terms in the series expansion (12); this reveals the dominant mode, with frequency 
equal to 2(jJt^ and the remaining smaller terms responsible for the anharmonicity in 
Figs. 1 and 2. Note that in the case of a pancake-like potential (D = 1), the well- 
defined parity of the transverse ground state requires to be orthogonal to states with 
odd parity. Therefore, the dominant mode in the decomposition (12) corresponds to the 
n = 2 term. 
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4. Adiabatic elimination of (p^ 

Due to the high anisotropy of the trapping potential, we expect the evolution of the 
longitudinal modes to take place on a much longer time scale than the transverse 
modes. In fact, our perturbation expansion is valid as long as this assumption holds, 
for according to Eq. (11), (pg = 0{e). In this section, we use this result to derive an 
adiabatic approximation for (p^, which we then use to bring Eq. (11) to a more tractable 
form. 

We start by integrating Eq. (10) by parts, which lets us write the following relation 
for (p^(t), n > 1: 

- -i^e-(^"-^o)(*-)/'^(^.„|cp^(.)r + ^.,|cp?(.)r)cp^(^^ \0{e). (19) 

I^n — -t^O ^ / 

Here the first term on the right-hand side arises from the straightforward integration of 
the rapidly oscillating exponential in Eq. (10), whereas the (neglected) remaining term 
is of higher order in e because it involves time derivatives of cpg and cpg*. We assume 
from now on, consistent with our perturbation theory, that these higher-order terms in 
Eq. (19) can be neglected. Under this approximation, we write 

<(t)^(p;^(t) + 4"(t), (20) 

where 

Kit) = -;|^ {9aamt)\' + ~9aMit)\') CP?W (21) 

corresponds to the adiabatic following of |(Pq pcpg and |(pQp(pQ, which is the dominant 
dynamical effect described by Eq. (19). The remaining lowest-order term, 

m = {<{o) - (p;^(o))e-(^"-^°)*/^ (22) 

is associated with the deviation of cp^(O), the projection of the initial condensate wave 
function onto (^^, from its adiabatic variant (p^(0). 

We can now use the adiabatic approximation (20) to eliminate (p^(t), n > 1, from 
Eq. (8), thus obtaining a more tractable evolution equation for cpg than Eq. (11): 

-26^Tt^,^(P^|(P^|2 (^^^^|cp^|2+^^J^^|2^ + ,2Jc.^ (23) 

where 



CXD 



n=l 

X 



(24) 
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n=l 



12 



(25) 



is determined solely by the properties of the transverse trap and characterizes the 
strength of the coupling of transverse and longitudinal directions. The quantity tj'^/Tt 
can be thought of as the relevant quantification of the transverse energy scale as far as 
the perturbation theory is concerned. 

After reorganizing the second-order terms in Eq. (23), we obtain the final form of 
the effective, coupled, two-mode evolution equation for the dominant longitudinal wave 
function (pQ, 

i/icp^ = e (^i^L +^aa^r I (poT +^a/3^r I cpoP) 



(26) 



where 



Here we introduce the functions 

oo 

= E(^-l^o)cPn(o)e-'(^''-^°)*/^ 

n=l 



vrit) = J2 



3\2 



MIL, 



-i(En-Eo)t/h 



(27) 

(28) 
(29) 
(30) 



Notice that Dt(0) = and Vrit) = irj^GT{t). 

The importance of the source term Iq in determining the dynamical behavior of (Pq 
essentially depends on the difference between (p"(0) and (p"(0). As shown by Eqs. (28) 
and (29), Iq is a rapidly varying function in comparison to the other terms in Eq. (26). 
In the case of a transverse harmonic confinement, we can write Vrit) in the closed form 



VT{t) 



Vt 

hjjJT 

Vt 



In 



ln( 1 - -e 



-2iujTt 



, D — 1 (pancake). 



D — 2 (cigar). 



(31) 



The absolute value and real and imaginary parts of Vt are plotted in Fig. 4 for D = 1 
and Fig. 5 for D = 2. 



Mean-field dynamics of two-mode BECs in highly anisotropic potentials 



13 




Figure 4. Plots of Dt, in units of {Jiwt / ^ as a function of time for a pancake-like 
geometry {D = 1) and transverse harmonic potential. Time is measured in units of 
uj^^. The dashed, black line represents the real part of Dt, the dotted, black line is 
the imaginary part, and the solid, purple line is the absolute value. 



0.15r 
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Figure 5. Plots of Dt, in units of {Jiwt / ^ as a function of time for a cigar- 
shaped trap {D = 2) with a transverse harmonic potential. Time is measured in units 
of uj^^ . The dashed, blue line represents the real part of Dt, the dotted, blue line is 
the imaginary part, and the solid, orange line is the absolute value. 
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^.7. Single-mode dynamics within the adiabatic elimination 0/ cp^ 

The dynamical equations for a single-mode BEC can be obtained trivially from Eqs. (20) 
and (26). By simply setting g^/s = and dropping the superfluous index a, we get 

(Pn = -^^^^glVoWo + In = Vn + L Tl > 1, (32) 

ihipo = e {Hl + gr]TW^\^ - 36^^Tt|(PoH (Po + (33) 

where 

in = ((Pn(0) - (p„(0))e-^(^»-^°)*/^ (34) 
/o = 2^|(pop($ - !•) + gipli^* - !•*), (35) 

and 

oo 

m = E(^-l^o)Vn(0)e-(^"-^°)*/^ (36) 

n=l 

m = -g\^)om'vMMt)- (37) 

Equations (32) and (33) correspond to the time-dependent versions of the stationary 
single-mode case [ ], but with the additional source terms and /q. In fact, in the 
stationary case, cp^ = (p^, and hence $ = l>. As a result, both and Iq vanish 
and the questions reduce to that of Ref. [ ]. Note that relative to a GP equation, the 
longitudinal equation (33) has an additional quintic term, which acts as an effective 
three-body, attractive interaction among the atoms regardless of the sign of g. This 
attractive interaction is mediated by the changes in the transverse wave function, as 
evidenced by the appearance of the (nonnegative) coupling parameter in the coupling 
strength Sg'^Tx- 

The coupling constants in Eq. (33) can be calculated explicitly for a transverse 
harmonic potential. It follows from Eqs. (16), (17), and (25) that 

2 

ln(8 - 4^3), D = 1 (pancake), 
^T=< , ^ (38) 
-^In-, i^ = 2 (cigar). 

As a result, the coupling constants for a cigar-shaped trap (rf = 1) are grjT = 2hujTa 
and 3g'^TT = 6/icJT<^^ ln(4/3), whereas for a quasi-2D pancake {d = 2), we obtain 
gr]T = 2V27rhujTpoa and Sg^Tr = 24.7rhwTpla^ ln(8 - 4^3). 

It is worth pointing out that such a self-focusing interaction has also been used 
to study the propagation of solitons in single-mode, attractive, quasi-lD condensates, 
trapped by an inflnitely long cylindrical harmonic potential [9, 10, 11]. In contrast to 
using a relative-state decomposition of the three-dimensional condensate wave function. 
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those studies use a Taylor expansion of a local, transverse chemical potential about 
the maximum density of the condensate to derive the additional quintic nonlinearity. 
This leads to an equation without the source term /q and with a coupling constant that 
differs from Sq^Tt (for = 1) by a factor of 4. This factor corresponds to the difference 
between using the maximum of the Gaussian ground-state probability distribution and 
the average ijt of the distribution over itself. 

5. Two- mode dynamics of a cigar-shaped ^'^Rb condensate 

In this section we apply the general results of the time-dependent relative-state 
decomposition derived above to analyze the two-mode dynamics of a cigar-shaped ^''Rb 
condensate consisting of atoms that can occupy the 5Si/2 |F = 1; = — 1) and 
|F = 2; Mf = +1) hyperfine states. These states are of particular interest to us 
because they can be used to implement a nonlinear Ramsey interferometry protocol 
with sensitivity that scales better than as discussed previously in Refs. [2, 4, 5, 6]. 

5.1. Two-mode evolution 

In typical experiments [12, 13, 14, 15, 16, 17, 18], the |F = 1; = — 1) = |1) state is 
trapped and cooled to the condensation point. Once the atoms in |1) have accumulated 
in the condensate ground state, a two-photon drive is used to couple the |1) state to 
the |F = 2; = +1) = |2) state. These two states have nearly identical magnetic 
moments, and hence feel essentially the same confining potentials. This strategy is used 
instead of cooling the two hyperfine states simultaneously to form a condensate in a 
superposition, because the lifetime of atoms in the |2) state in a trap is much shorter 
than the lifetime of atoms in the |1) state. Thus, we assume from now on that the optical 
pulse suddenly creates the superposition state '0(r)(|l) + \2))/^/2 for each atom.:j: This 
procedure creates the two-mode condensate with both modes being occupied by the 
same number of atoms and having the same initial wave function ip{r). In this case, the 
inter-mode coupling is symmetric, i.e., gi2 = ^21 = guiN — l)/2, and we can write the 
longitudinal two-mode equations (26) as 

ih(pl = e{HL + guVrlvll^ + 9i2Vt\vI\^) cPo 

- e^TT{3gl,\ipl\^ + (4^n + 2^12)^12 Icp^HcPor + ^12(^12 + 2^22)|cp^|') 

+ e'/o, (39) 

J This approximation holds as long as the Rabi-flopping time scale for the atomic transition is smaller 
than the time scale for the transverse BEC dynamics (the longitudinal BEC dynamics occurs on an 
even longer time scale). For the ^^Rb BECs we consider here, the Rabi-flopping time scale is typically 
less than a millisecond, so these conditions hold for the trapping potentials we consider here. 
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ihdpl = e{HL + g22r}TWl? + 9i2'nTWl?) ^>l 

- e^TT(3^22l<Pol^ + (4^22 + 2^12)^12 1 (PoHtPoT + ^12(^12 + 2^ii)|(Po|^)(Po 
+ e'/o , (40) 

which, of course, are written under the assumption that the longitudinal functions <^\{t) 
and (p^(t) can be approximated by Eq. (20) and, hence, are given by 

<PnW = (^iikowr + -mwim <PoW + (41) 

— ^0 

<Pl{t) = {~92Mit)\' + ~9iMm <pl{t) + ilit). (42) 

Before we write down the expressions for the source terms in the equations above, we 
use below the results from the relative-state decomposition of the ground-state wave 
function, derived in [ ], to bring /q, /q, and to a simpler form. 

In our mean-field model, we suppose that all atoms initially occupy the same 
single-particle wave function r), which is the solution of the time-independent, 
three-dimensional GP equation for the state |1). For this particular choice of the initial 
state, we showed in [ ] that, in the perturbative regime <C A/r, the condensate 
wave function can be approximated, to first order in perturbation theory, by the time- 
independent relative-state decomposition 

oo 

ij{p,r) = Up)Mr) + eJ2Up)Mr) = i^\t = 0) = ij'{t = 0), (43) 

n=l 

where the dominant longitudinal wave function (foi'f^) is the solution to 

(Al -Hl- gVT^l + 3efTT^^,)^o = 0, (44) 
with g = {N — l)gu^ and 

^no = -gv^loP^T = <p\it = 0) = <pl{t = 0). (45) 

Here (foo is the zero-order longitudinal wave function given by the reduced-dimension, 
single-mode GP equation 

(Ml - Hl- gr]T^m)^m = 0. (46) 
Thus, using Eq. (45), we can write Eq. (20) as 

fa{-t\ — (7, 7, 7, ^V^3 (^olCn) -i(En-Eo)t/h ( 



and the functions (28) and (29) as 

^\t) = ^\t) = -g^l,VT{t), (48) 
^-P{t) = - {~g^^ + ~g^p) <floMt)- (49) 
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In Eq. (49), we used that (Po(0) = 9:^00 + 0{e) and kept only the lowest-order terms, 
considering that the source terms are of quadratic order in e. From Eqs. (48) and (49), 
it follows trivially that 

$1 _ = _ g^^ _ g,,)^l^VT{t), (50) 

$2 _ 1,21 _ _ _ g^^)^l^v^(t). (51) 

For our particular choice of the hyperfine levels of ^''Rb, the 5-wave scattering 
lengths for the processes |1)|1) ^ |1)|1), |1)|2) ^ |1)|2), and |2)|2) ^ |2)|2), 
respectively, are an = 100.40ao, = 97.66ao, and = 95.00ao [1 with qq being the 
Bohr radius. Thus, the inter-mode coupling is given approximately by the arithmetic 
mean of the intra-mode coefficients, i.e., gi2 ^ {gu + fl'22)/2, which equivalently means 
that 

^11 - ^12 - (^11 - ^22)/2 = 7i- (52) 

Therefore, it follows that 

9-911- 912 = 1(911 - 9i2){N - 1) « 271 (iV - 1), (53) 
9 - 922 - 912 = ^{2911 - 912 - 922){N -l)^ 371 (A^ - 1). (54) 

We can now put all these results together with Eqs. (27) and (47) to write all the 
source terms as 



/o^(t) = -27i(iV - 1)^3^1 (2^n|cPoWr + ^i2|<PoWI') VT(t) +^ii(<Po(t))'t);.(t) 

+ \g,2^l{t) [^>l{t)v*r,{t) + ^>l*{t)vT{t)\ }, (55) 

ilit) = -37l(iV - l)iplo{ (2^22|cp^(t)r +^12|(pJ(t)r)l^T(t) +^22((P^(t))'-UTW 

+ l~9i2<pl{t) Wl{t)vUt) + <Po{t)vT{t)] }, (56) 
m = -27i(iV - l)^^o#^e-(^»-^o)*/^ (57) 

P^it) = -37i(iV - l)^3^ity^e-(^"-^o)*/^ (58) 

— ^0 

Because 71 <C 5^22, 5'i2, S'li, however, we expect the contribution from all the source terms 
to the two-mode evolution to be negligible. In other words, since the initial longitudinal 
function (fn{0) is not much different from cp^, the evolution of cp^ should not differ much 
from the adiabatically tracking solutions. 



5.2. Numerical simulations 

To study the two-mode evolution described by our perturbative model, we numerically 
integrate Eqs. (39) and (40) and compare it against the exact numerical results for the 
evolution of the two-mode BEC, given by the two-mode, coupled, three-dimensional 
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Gross-Pit aevskii equations (1) [ ]. We restrict our numerical analysis to the same 
scenario of the three-dimensional simulations presented in Ref. [ ]. That is, we consider 
the two-mode condensate to be trapped by cigar-shaped potentials of the form 

V{p,z) = ]^{mujlp' + kz'^), (59) 

with g = 2, 4, and 10. These three potentials allow us to explore how the results vary 
due to the inhomogeneity of the trapping potentials. Notice that the limit q ^ oo 
recovers the case of a homogeneous, hard- walled trap. 

We set the transverse frequency to 350 Hz and the longitudinal frequency to 3.5 
Hz for the harmonic case (g = 2), with the result that Mt — 14 000 atoms. We obtain 
this estimate of Nt from the definition given in [ ] . To compare the simulations for the 
different power-law potentials, we choose the stiffness parameter k so that Mt remains 
the same for the two other values of g; thus all the traps have the same one-dimensional 
regime of atom numbers. For such choice of parameters, we find po — 0.6 /xm and the 
aspect ratio of the bare traps (po^^o) to be approximately equal to 1:10, 1:24, and 1:57, 
respectively, for g = 2, 4, 10. 

5.3. Differential phase dynamics: nonlinear Ramsey interferometry 

We are now in position to revisit the numerical simulations of a nonlinear Ramsey 
interferometer presented in Ref. [ ] and use the results of the time-dependent relative- 
state decomposition to analyze the relative-phase dynamics of the two-mode condensate. 

As in typical Ramsey interferometry schemes, the protocol runs as follows. The 
atoms are first condensed to the state ?/;(r)|l), and a fast optical pulse suddenly creates 
the superposition state '0(r)(|l) + \ 2))/^/2 for each atom. The atoms are then allowed to 
evolve freely for a time t, which brings the atomic state to ['^'^(r, t)|l)+'0^(r, t)|2)]/^/2. A 
second transition between the hyperfine levels is then used to transform any differential 
phases between the two modes into population information that is finally detected. For 
this second transition, we choose a 7r/2 rotation about the Bloch x axis, changing the 
atomic state to 

_(^i_,^2)|i^_^(^i+^^2)|2). (60) 

Thus any differential phases accumulated by the mean-field solutions will give rise to 
an interference fringe pattern in the detection probabilities for each hyperfine level, 

p,^, = \[lTlm{m^'))\, (61) 
due to the overlap of the two spatial wave functions, 

(t/;^!^/;^) = j d^rij^'^\ (62) 

Note that the imaginary part of the overlap (62) is responsible for the fringe pattern in 
this interferometry scheme. 
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5.4' Time scales 

The relevant time scale for the relative-phase dynamics of the two-mode condensate can 
easily be estimated if one completely neglects the spatial evolution of the condensate. 
Under this approximation, the two-mode evolution becomes trivial. Due to the difference 
in scattering lengths, the only effect of the evolution is to introduce a differential phase 
shift between the states |1) and |2). This relative phase simply corresponds to the 
difference between the average scattering energies of the condensate, which implies that 
the detection probabilities (61) oscillate as 

Pi^2 = \[lTM^Nt)\, (63) 

where 

fijv = (iV - l)vii/h (64) 
is the fringe frequency. Here the quantity 



r/ = y rfV|V^(r)|^ (65) 

is a measure of the inverse volume occupied by the ground-state wave function. This 
fringe pattern allows one to estimate the coupling constant 71 with an uncertainty given 
by 

\d{Jy)/d^,\ ^/NNr^' 

whose scaling with the atom number is better than the optimal 1/A^ scaling of linear 
interferometry [5]. 

Within the reduced-dimension approximation of Eq. (3), the dominant contribution 
to the fringe frequency is 

^rda = Nr]TVL7i/fi, (67) 

where 

Vl= f d^r\^oo{r)\\ (68) 



This frequency sets the characteristic time scale of the short-time relative evolution of 
the longitudinal wave functions cpj and (pg- 

In Fig. 6, we plot the values of the Ramsey fringe frequencies and f^rda, as 
a function of the atom number and the three different values of q. We calculate 
using the numerical evaluation of rj. Not surprisingly, f^rda gives an upper bound on 
the numerical values oi Q^. The important point here is not the difference between 

and f^rda, but rather that both are much smaller than the transverse frequency 
ujt^ as is to be expected for the highly anisotropic potentials used in our simulations. 
This result indicates that the longitudinal frequencies do satisfy the adiabatic condition 
Qn^ f^rda <^ ^T- We therefore expect that the longitudinal two-mode equations (26) can 
be used to model the dynamics of the BEC interferometry protocol. 
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Figure 6. Idealized fringe frequency Vt^ = Nrjji/h, in units of the bare transverse 
frequency cjt, as a function of the number of condensed atoms. This fringe frequency, 
calculated using the numerical results for t^, is plotted as circles (blue) for (7 = 2, squares 
(black) for ^ = 4, and triangles (red) for q = 10. The dashed lines correspond to the 
reduced-dimension estimate, flrda = NtjttjlJi/^^ which gives an upper-bound to the 
fringe frequency Qn calculated using the numerical results for rj. Both frequencies are 
much smaller than the transverse trap frequency ujt- 



5.4-1' Ramsey fringes. According to the relative-state decomposition (6), it is easy 
to see that, at the order we are working, the dominant longitudinal functions (pg and 
(Po carry all the information about the relative phase between the three-dimensional 
condensate wave functions t/^^ and t/^^, since 

m^') = {^>lWl)+0{e'). (69) 

By using the numerical solutions of the time-dependent, coupled, 3D GP equations, we 
calculate the spatial overlap and by using the numerical solutions of Eqs. (39) 

and (40), we calculate the overlap (cpol^Po)? both as a function of time for various atom 
numbers. In addition, we calculate the spatial overlap between the numerical solutions 
of the quasi- ID, time-dependent GP equation (the reduced-dimension approximation) 
to use as a benchmark for our (relative-state-decomposition perturbation) model [ ]. 
We stress that these spatial overlaps and the resulting Ramsey fringes, including 
those coming from the reduced-dimension approximation, are more than the simple 
estimates (64) and (67), since the numerical integrations include effects of position- 
dependent phase shifts and reduction in fringe visibility (see Sec. 5.4.3) not captured by 
the simple estimates of fringe frequency. 

In Fig. 7, we compare the Ramsey fringes obtained from the imaginary part of 
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time (ms) 

Figure 7. Ramsey fringes for a cigar-shaped ^^Rb BEC of 500 atoms. The points 
correspond to the results of the numerical integration of the time-dependent, coupled, 
two-mode, three-dimensional GP equations (4.22) for the different trapping potentials 
(4.13): circles (blue) correspond io q = 2^ squares (black) to = 4, and triangles (red) 
to g = 10. The corresponding solid lines are the numerical results of the integration of 
Eqs. (39) and (40), given by the relative-state-decomposition perturbation theory. The 
dashed lines represent the results coming from the numerical integration of the time- 
dependent, two-mode, quasi-lD GP equation (the reduced-dimension approximation). 
For all three trap geometries, the agreement between our perturbation theory and the 
exact 3D numerics is remarkably good during the whole Is of integration time. The 
reduced-dimension approximation, on the other hand, only describes earlier stages of 
the evolution. As q increases, both models achieve a better performance. Note that 
the reduction in the fringe visibility is small, being almost absent for q = 10. 

these three spatial overlaps for a condensate of 500 atoms. For all three longitudinal 
potentials, our perturbation theory reproduces the exact 3D numerics remarkably well 
during the entire integration time interval. The reduced-dimension approximation, on 
the other hand, only holds in the earliest stages of the evolution, for about a quarter of a 
fringe period. It is no surprise that as the traps get harder and more homogeneous, both 
models achieve a better performance. Note that in this case of rather small nonlinearities 
{N/Mt ^ 3.5 X 10~^), the reduction in the fringe visibility due to position-dependent 
differential phases is quite small, being almost absent for q = 10. This effect becomes 
obvious for stronger couplings, however, as illustrated in Fig. 8 for a condensate of 1 000 
atoms. For this many atoms and same time window of 1 s, the agreement between our 
perturbation model and the 3D solution is still quite good, especially for the more 
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Figure 8. Ramsey fringes for a cigar-shaped ^^Rb BEC of 1 000 atoms (labeling 
convention as in Fig. 7). The points represent the numerical results of the integration 
of the time-dependent, coupled, two- mode 3D GP equations for the different trapping 
potentials (59), whereas the lines are the relative-state (solid) and quasi-lD (dashed) 
predictions. The agreement between our model and the 3D solution is quite good, 
specially for the g = 4 and (7 = 10 cases. Both fringe frequency and amplitude are better 
predicted by our model than by the reduced-dimension (quasi- ID) approximation. For 
q = 2^ small deviations from the expected signal accumulate over time, resulting in a 
phase mismatch around t = 500 ms. 

homogeneous trapping potentials g = 4 and q = 10. Both fringe frequency and 
amplitude are well predicted by our model, showing a substantial improvement over 
the reduced-dimension approximation. For a harmonic trap, small deviations from the 
expected signal accumulate over time, resulting in an evident phase mismatch around 
t = 500 ms. 

For a condensate of 2 000 atoms, our model starts to break down, as shown in Fig. 9. 
For all g's, the reduction in fringe visibility is still well captured by our perturbative 
equations, with agreement with the 3D numerics getting better as q increases, but the 
predicted fringe frequency is a bit too small. This effect becomes more obvious for a 
larger condensate with 5 000 atoms. As clearly seen in Fig. 10, our prediction of the 
fringe frequency is too small and thus shows the complete breakdown of our perturbation 
model, which fails to describe the correct relative-phase accumulation. 
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Figure 9. Ramsey fringes for a cigar-shaped ^^Rb BEC of 2 000 atoms trapped by the 
potentials (59) with (a) q = 2^ {h) q = 4, (c) q = 10. The points represent the numerical 
results of the integration of the time-dependent, coupled, two-mode 3D GP equations. 
Solid lines are the corresponding relative-state predictions, whereas the dashed lines 
are results coming from the numerical solution of the quasi- ID GP equation. The 
reduction in fringe visibility is still well captured by our model, with agreement with 
the 3D numerics getting better as q increases. The predicted fringe frequency, however, 
is a bit too small, which points to the breakdown of the perturbation theory, yet more 
evident for 5 000 atoms in Fig. 10. 



5.4-2. Source terms. As a brief remark, we point out that in all the results presented 
in this section, the solutions to Eqs. (39) and (40) take into account the source terms 
that arise due to the difference cp^(O) — (p^(0). We do expect that the contribution 
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Figure 10. Ramsey fringes for a cigar-shaped ^^Rb BEC of 5 000 atoms trapped 
by the potentials (59) with (a) g = 2, (b) g = 4, (c) q = 10. The points represent 
the numerical results of the integration of the time-dependent, coupled, two-mode 3D 
GP equations. Solid lines are the corresponding relative- state predictions, whereas 
the dashed lines are results coming from the numerical solution of the quasi- ID GP 
equation. The perturbation theory breaks down: Eqs. (39) and (40) can no longer 
correctly describe the relative-phase accumulation by the condensate mean fields. 



of the source terms is quite small for the initial condition we use. To verify this, we 
plot in Figs. 11 and 12 the predicted Ramsey fringes for 1000 and 5 000 atoms using 
the solutions of Eqs. (39) and (40) with and without the contribution from the source 
terms. The difference in the computed overlaps is indeed quite small, being of order 
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Figure 11. Ramsey fringes for 1000 atoms using the solutions of Eqs. (39) and (40) 
with (sohd hue) and without (dash-dotted Une) the contribution from the source terms 
for g = 2,4, 10 (blue, black, red). The two lines lie almost on top of one another. The 
inset shows a closeup for the q = 2 case around t = 465 ms, showing a difference of 
the order of 10~^ between the two fringes. 

10~^ for the case of 1 000 atoms and 10~^ for 5 000 atoms, as shown in the insets of 
Figs. 11 and 12 for the harmonic trap. 

5.4-3. Thomas- Fermi estimates. For the traps and atom numbers that we are 
considering, it is legitimate to ignore the kinetic-energy term in Eqs. (39) and (40) 
and work in the longitudinal Thomas-Fermi approximation. Within this regime, the 
probability densities do not change with time, i.e., 

|(p^(z,t)P = |(po(z,0)|2, (70) 

with |(po(2;, 0)p given by the Thomas- Fermi solution of Eq. (44), 

(Po(^,0) =— + e [— =Qo{z). (71) 

Under this approximation the two-mode evolution can be described by a simple 
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Figure 12. Ramsey fringes for 5 000 atoms using the solutions of Eqs. (39) and (40) 
with (sohd hue) and without (dash-dotted Une) the contribution from the source terms 
for g = 2,4, 10 (blue, black, red). The two lines lie almost on top of one another. The 
inset shows a closeup for the q = 2 case around t = 415 ms, showing a difference of 
the order of 10~^ between the two fringes. 



accumulation of a phase that depends on the local atomic linear density, given by 

it (\ 



n\2 



-kz"^ + r)TQo{z){gaa + 



This yields an overlap 

where the position-dependent differential phase shift is 

oO{z) = 1 - e — Qo{z){3gu + ^922 + ^9x2, 

n \ vt 



■ (72) 



(73) 



(74) 



Note that we can easily retrieve the Thomas-Fermi prediction of the quasi- ID 
approximation by taking the limit e ^ 0. This gives the following approximation for 
the overlap (62) 

iij'lij') ^ {4>'\(l>') = [ dz qo{z)e-'''^^\ (75) 
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Figure 13. Same as Fig. 8 (1 000 atoms), but here the hnes correspond to the Thomas- 
Fermi estimate (73) for q = 2,4,10 (blue, black, red). Equation (73) provides a very 
good account of the exact 3D fringe signal (dots), with better agreement for short times 
and for harder traps. For longer times, however, the analytical Thomas- Fermi model 
breaks down, which suggests that the probability densities can no longer be regarded 
as constants in time. 



where the relative phase is 



5e{z) = - l)^^go(^)^i ^ (76) 
h 



with qo{z) being given by 

fiL - kzy2 

qo[z) = z . 77 

9Vt 

Here /x^, the zero-order longitudinal part of the chemical potential, is determined from 
the normalization condition for qo{z)^ which gives 

Note that, as far as the spatial overlap (73) is concerned, not only the transverse, 
but also the longitudinal correction play a role in rectifying the reduced-dimension 
approximation (75). Comparing the two overlaps, one can interpret the longitudinal 
correction as coming from the distribution Qq{z)^ whereas the transverse correction 
comes from a rectified tjt^ which is given by 



rir(^) =Vt- eTTQo{z){3gu + 3^22 + 2^12). 



(79) 
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Figure 14. Same as Fig. 9 (2 000 atoms), but here the hnes correspond to the Thomas- 
Fermi estimate (73) for g = 2, 4, 10 (blue, black, red). The agreement between the exact 
3D overlap and Eq. (73) is still good for the harder traps (7 = 4 and q = 10. For q = 10, 
in particular, the mode provides the correct prediction for the entire integration time. 
For ^ = 2, however, the model quickly deviates from the expected 3D numerical signal 
(dots), due to changes in the probability densities. 

Note that r[T{z) is a function of the local axial distribution Qo{z) and, hence, takes into 
account inhomogeneities in the trapping confinement. 

In Fig. 13, we show the imaginary part of the overlap (73) for 1 000 atoms and 
compare it with the exact 3D overlap and the reduced-dimension prediction given by 
Eq. (75). The Thomas- Fermi estimate (73) provides a very good account of the exact 3D 
fringe signal, being far superior to the quasi- ID model. As expected, the agreement with 
the 3D numerics gets better for higher q and shorter times. For longer times, however, 
the assumption of a constant probability density is no longer valid, and the Thomas- 
Fermi approximation breaks down. Figure 14 shows the same comparison for 2 000 
atoms, but without the ineffective quasi- ID results. For q = 2 and q = 4:^ deviations 
from the exact fringe pattern become significant after a quarter of a fringe period. For 
q = 10, however, our analytical prediction using the Thomas-Fermi approximation is 
remarkably accurate. 

5.5. Average densities 

To examine further the extent to which our perturbative approach correctly describes the 
two-mode evolution, we study below the dynamics of the average mean-field densities for 
various atom numbers. In particular, we can use the perturbation theory to describe the 
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density evolution up to a point where the spatial separation becomes clearly apparent 
(see Fig. 15). 

According to the perturbative expansion (6), the average single-mode density per 
atom, to first order in e, is given by 

ut) = J d\\r{t)\' (80) 

oo « 

= VTVl{t) + 2eJ2{Ueo) / c?^|(p^(t)P((prW<W + c.c.), (81) 

n=l ^ 



where 



V2it) = J dz\ip^{t)\\ (82) 

Within the adiabatic approximation (20) to (p^(t), we can split Eq. (81) into two 
terms, a dominant term and a second term that comes from the source term in Eq. (20): 

Ve.{t)=Ut) + f)f{t). (83) 
The dominant term, which is given by 

oo « 

Ut) = VTV2{t) + 2eJ2{^^\Q / dz\<p^{t)W<prm^{t) +C.C.), (84) 

n=l ^ 

includes the contributions that arise from (Pq and cp^. By using Eq. (21) and introducing 
the quantities 

Vteit)- J dz\ip^{t)\', (85) 

vr\t)= I d^i(p-(t)n<p?(t)r, (86) 

we can write the dominant term in the more compact form 

f)a{t) = VTVlit) - 4€Tr {gaaVUt) + ~9.pVT^ {t)) . (87) 

The second term in Eq. (83), 

oo « 

fin) = '2^Y.^in\ii) J dz\<p^it)\'{<prm{t)+c.c.), 



n=l 



comes from the source term I^. Together with Eqs. (57) and (58), this gives 

r7}(t) = -467i(iV-l) J dzcpl,\^>l{t)\'l^^l*{t)vT{t) + c.c}j, (89) 



fij{t) = -6ej,{N-l) j dz^lM{t)\''[ipl\t)vT{t) + c.c. 



(90) 
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Figure 15. Density evolution for a BEC of 1 000 atoms and q = 2. The subplots 
on the left show the time evolution of whereas the time evolution of 

t)p is illustrated by the figures on the right. Initially, both modes have the 
same density profile as seen in subplots (a) and (d). The position-dependent phases 
acquired by each mode drive differences between the atomic densities associated with 
the two hyperfine levels, as seen in subplots (b) and (e) for 320 ms of evolution. Atoms 
in mode 1 are pushed towards the edges of the trap due to the stronger intra-mode 
repulsion (^n > gi2 > ^22), while atoms in mode 2 get compressed in the center of 
trap by the atoms in mode 1 because of the inter- mode repulsion {gi2 > 922)- This 
leads to spatial separation of the two modes, as shown in subplots (c) and (f) at 700 ms 
of the evolution. Due to the trap confinement, this dynamics is eventually reversed. 
The dynamical relative- state perturbation theory allows us to investigate the density 
evolution long enough to capture all of these effects. 
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A similar set of steps, using Eqs. (6) and (20), splits the density overlap, 

Vi2{t) = I d'r\^\t)f\^\t)\^ (91) 

into two terms, 

Vi2{t) = m2{t)+fi}^{t), (92) 

where 

Vuit) = VTVLit) - 2eTT[{~9n + ^12)^7l + + ^i2)r?f '(t)], (93) 

= I dz\<pl{t)WM', (94) 

and 

77f(t) = -3e7i(iV-l) j dz^lM{t)\'(^<pl*{t)vT{t) + c.c}j 

-2eJ^{N-l) J dz^l,\<pl{t)\'(^^l*it)vT{t) + c.c}j. (95) 

Our first step in studying the time evolution of the mean-field densities is simply to 
calculate ^2(^)5 and 7712 (i), as instructed by Eqs. (80) and (91), using the numerical 

solutions of the time-dependent 3D GP equation (5). We compare these to the relative- 
state predictions (83) and (92). Again, we also calculate the corresponding predictions of 
the reduced-dimension approximation to use as a benchmark to our model. Figures 16, 
17, and 18 show these three quantities in harmonic trap units [ ] for g = 2,4 and 10. 
In each figure, we consider the cases for 500, 1000, 2 000, and 5 000 atoms. We note 
that the 3D numerical solutions, the relative-state model, and the reduced-dimension 
approximation all use different initial conditions. 

In all cases, we can see that the 3D mean densities are poorly described by 
the reduced-dimension approximation. Our perturbative model, however, manages to 
reproduce the exact mean densities remarkably well for 500 and 1 000 atoms. Small 
deviations appear for 2 000 atoms; for 5 000 atoms, the model gives entirely incorrect 
predictions for the density evolution. The case q = 10 and = 500 is an exception 
to the previous statement. For this particular case, numerical instabilities affected 
the calculation of the initial longitudinal wave functions cpo(O) and (Po(0), which 
consequently resulted in the inaccurate density dynamics seen in Fig. 18(a). Note, 
however, that this numerical error does not seem to affect the good agreement for the 
relative phases shown in Fig. 7. 

The two-mode density evolution shown in Figs. 16-18 can roughly be described as 
follows. Initially all densities are the same, i.e., r]i{0) = 772(0) = 7712(0); hence, the earlier 
stage of the dynamics is essentially dictated by the different scattering couplings. Thus, 
atoms in mode 1 are initially pushed towards the edges of the trap due to the stronger 
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Figure 16. Time evolution of the average densities, rji (orange) and 772 (green), and 
of the density overlap, 7712 (purple), for = 2 and (a) N = 500, (b) N = 1000, 
(c) N = 2 000, and (d) N = 5 000. The points represent the numerical results of 
the integration of the time-dependent, coupled, two-mode 3D GP equations. Solid 
lines are the respective relative-state predictions, whereas the dashed lines are results 
coming from the numerical solution of the quasi- ID GP equation (reduced-dimension 
approximation). In all cases, the 3D mean densities are poorly described by the 
reduced-dimension approximation. Our perturbative model, on the other hand, 
performs extremely well, except for = 5 000, where the model breaks down. 



intra-mode repulsion {gu > gi2 > 5^22) , while atoms in mode 2 get compressed in the 
center of trap by the atoms in mode 1 because of the inter-mode repulsion {g^ > 5^22) • 
Therefore, the wave function for mode 1 gets wider than the wave function of the second 
mode. In other words, the mean atomic density rji decreases while r]2 increases, and as 
the two modes are driven apart, the density overlap rju drops. Of course, the atomic 
repulsion is counterbalanced by the trapping potential, so at some point this dynamics 
is reversed. 

The complex nonlinear evolution in Figs. 16-18 makes it challenging to analyze 
further the dynamical behavior of the mean-field densities. To gain a deeper 
understanding of the results, we examine the renormalized single-mode mean densities. 
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Figure 17. Time evolution of the average densities, rji (orange) and 772 (green), and 
of the density overlap, 7712 (purple), ioi q = A and (a) N = 500, (b) N = 1000, 
(c) N = 2 000, and (d) N = 5 000. The points represent the numerical results of 
the integration of the time-dependent, coupled, two-mode 3D GP equations. Solid 
lines are the respective relative-state predictions, whereas the dashed lines are results 
coming from the numerical solution of the quasi- ID GP equation. Our perturbative 
model performs extremely in all cases, except for N = 5 000, when the model breaks 
down. Note that for 2 000 atoms, small deviations appear. The reduced-dimension 
approximation fails to describe the 3D mean densities in all cases. 



r]i and 772, which we define according to 

= ^, (96) 
V 

where 

V = — ^ — (97) 

is the arithmetic mean of the two average densities. 

Figures 19, 20, and 21 show the time evolution of 7)1 and 7)2 for 9 = 2,4, and 10. In 
each figure, we plot the cases = 500, 1000, 2 000, and 5 000. The time evolution of 
771 and 7/2 is, in fact, much simpler than of rji and 772, and, therefore, easier to interpret. 
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Figure 18. Time evolution of the average densities, rji (orange) and 772 (green), and 
of the density overlap, 7712 (purple), for g = 10 and (a) N = 500, (b) N = 1000, 
(c) N = 2 000, and (d) N = 5 000. The points represent the numerical results of 
the integration of the time-dependent, coupled, two-mode 3D GP equations. Solid 
lines are the respective relative-state predictions, whereas the dashed lines are results 
coming from the numerical solution of the quasi- ID GP equation. The inaccurate 
density dynamics seen in subplot (a) is due to numerical instabilities that affected the 
calculation of the initial longitudinal wave functions cpo(O) and cpo(O). This numerical 
error, however, does not seem to affect the good agreement for the relative phases 
shown in Fig. 7. 



Among other things, using the renormalized average densities reduces the effect of the 
differences in initial conditions seen in Figs. 16-18. For instance, it is easy to see in 
Figs. 19-21 both the symmetry in the relative dynamics of the two modes and the 
suppression of the spreading of the condensate wave function with increasing q. Indeed, 
note that the harder the trap, the closer f]i and 7)2 get to 1. 

More importantly, the comparison of the dynamics of rj^ with of fj^ reveals an 
interesting fact about the dynamics of our model. The results for the renormalized 
average densities are clearly better described by our model than the evolution of rji and 
r]2. The cases of a condensate with 5 000 atoms, shown in subplot (d) of Figs. 16-18, is 
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Figure 19. Time evolution of the renormalized average densities, rji (orange) and 772 
(green), for (7 = 2 and (a) N = 500, (b) N = 1000, (c) N = 2 000, and (d) N = 5 000. 
The points represent the numerical results of the integration of the time-dependent, 
coupled, two-mode 3D GP equations. Solid lines are the respective relative-state 
predictions, whereas the dashed lines are results coming from the numerical solution 
of the quasi- ID GP equation. 

a good example of what we mean. In these three cases, our model failed to predict the 
correct density evolution. When we analyze the renormalized densities in subplot (d) of 
Figs. 19-21, however, this is clearly no longer the case. The agreement between the exact 
3D results and our perturbation model is overall reasonably good, and it gets better for 
shorter times and harder traps, as is to be expected. Surprisingly, the same can also be 
said for the reduced-dimension predictions, for there is a dramatic improvement in the 
agreement between the 3D exact dynamics and the quasi-lD model. In fact, the reduced- 
dimension approximation failed in all the cases shown in Figs. 16-18. Figures 19-21, 
on the other hand, show that the evolution of f]i and 7)2 can be well described by the 
quasi- ID GP equation, specially for smaller atom numbers. These results are indeed 
very revealing, for they clearly indicate that Eqs. (39) and (40) predict, to some extent, 
the correct relative density dynamics in spite of the use of an imprecise initial condition. 




Figure 20. Time evolution of the renormalized average densities, r]i (orange) and 772 
(green), for (7 = 4 and (a) N = 500, (b) N = 1000, (c) N = 2 000, and (d) N = 5 000. 
The points represent the numerical results of the integration of the time-dependent, 
coupled, two-mode 3D GP equations. Solid lines are the respective relative-state 
predictions, whereas the dashed lines are results coming from the numerical solution 
of the quasi- ID GP equation. 



which is ultimately related to underestimating the transverse profile of the condensate. 

Finally, we introduce the renormalized density overlap, J^^, as the ratio of 7712 to 
the geometric mean of the single-mode average densities, i.e., 

= (98) 

This quantity is equal to zero if there is no overlap between the mean-field distributions 
and is equal to one when they completely overlap, thus working analogously to the 
fidelity of quantum states. Figure 22 shows the time evolution of Eq. (98) for = 500, 
1000, 2 000, and 5 000 and the three different traps g = 2,4 and 10. This figure nicely 
shows the time scale on which the condensate density is constant, which increases as the 
trap gets harder, as expected. During this time interval, the two-mode dynamics simply 
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Figure 21. Time evolution of the renormalized average densities, 771 (orange) and 772 
(green), for (7 = 10 and (a) N = 500, (h) N = 1 000, (c) N = 2 000, and {d) N = b 000. 
The points represent the numerical results of the integration of the time-dependent, 
coupled, two-mode 3D GP equations. Solid lines are the respective relative-state 
predictions, whereas the dashed lines are results coming from the numerical solution 
of the quasi- ID GP equation. 



correspond to the accumulation of a position-dependent phase-shift. In addition, the 
agreement between our model and the 3D numerics is very good for the entire integration 
time, except for 5 000 atoms, where we see the model breaking down around 200 ms. 
This result enforces our previous conclusion that Eqs. (39) and (40) predict to some 
extent the correct relative density dynamics in spite of the use of an incorrect initial 
condition. 

6. Optimal reduced-dimension equations for the two-mode evolution 

Having successfully obtained a perturbative description of the time-dependent mean 
fields of the two-mode BEC, we turn now to the question of how to decouple the tightly 
confined dimensions from the mean-field evolution. 
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Figure 22. Time evolution of J"^ = 7?i2/vW?2 for (a) N = 500, (b) N = 1000, 
(c) N = 2 000, and (d) N = 5000 and the three different traps g = 2,4, and 10 (blue, 
black, and red). This quantity is zero when there is no spatial overlap between the 
mean- field densities and is equal to one when they completely overlap. The points 
represent the numerical results of the integration of the time-dependent, coupled, 
two-mode 3D GP equations. Solid lines are the respective relative-state predictions, 
whereas the dashed lines are results coming from the numerical solution of the quasi- ID 
GP equation. The agreement between our model and the 3D numerics is very good for 
the entire integration time, except for 5 000 atoms, where we see the modeling breaking 
down around 200 ms. 

For a single-mode BEC, we showed in [^] that, in the stationary case, such 
a decoupling can be formally achieved by performing a perturbative Schmidt 
decomposition of the ground-state condensate wave function, r), between the 
transverse and the longitudinal directions. In the Schmidt perturbation theory, we 
seek a solution to the time-independent, single-mode GP equation in the form 

r) = V>^oXo{p)Mr) + e^,Xi{p)Mr) + 0{e'), (99) 

where Xn{p) and (/>n('^) are the orthonormal Schmidt functions in the transverse and 
longitudinal directions and the vXi's are the (nonnegative) Schmidt coefficients (the 
squares are the eigenvalues of the marginal transverse and longitudinal density 
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matrices) [ ] . Note that we treat the Schmidt decomposition formaUy as a power-series 
expansion in the perturbation parameter e. Thus, the dominant term in Eq. (99) gives 
the best product approximation to the exact condensate wave function [2z], whereas the 
remaining terms introduce nonseparable corrections and, hence, describe entanglement 
between the transverse and longitudinal degrees of freedom. Because in the perturbative 
regime the spatial entanglement is small [ ] , one can neglect the nonseparable terms and 
describe the condensate by a product wave function. In this way, the transverse and 
longitudinal directions are decoupled and the dominant longitudinal Schmidt function 
gives the optimal reduced-dimension description for the condensate mean fields. Below 
we extend this formalism to the time-dependent mean fields of the two-component BEC, 
by showing that the instantaneous perturbative Schmidt decomposition can be obtained 
from the relative-state decomposition (6). 

In the single-mode case, the mean field (99) has the form of a two-qubit state to 
first order in e. In the two-mode case, however, the single-particle condensate state has 
the following form 



which corresponds to a tripartite state, for it has an additional qubit degree of freedom, 
namely, the two possible hyperfine configurations. Because the Schmidt decomposition 
only applies to bipartite systems, the decomposition of tripartite states is not unique 
and depends on the specific way one chooses to partition the system. 

From all possible bipartitions, the one that is physically relevant for understanding 
the reduced-dimension structure of the mean field of a highly anisotropic condensate 
is the one that separates the transverse degrees of freedom from the longitudinal and 
hyperfine (internal) degrees of freedom. This partition relies on the physical difference 
between the high energy scale of the tightly confined dimensions and the much lower 
energy scale of the longitudinal and hyperfine degrees of freedom. 

The procedure we describe next is very similar to the one discussed in [ ] . We start 
by pointing out that the two-mode, single-particle state (100) can be written in the form 
of the spinor 



which represents a vector in the two-dimensional hyperfine space, whose components 
are the condensate wave functions 



p,r, t)|2)] IV2., 



(100) 




(101) 




(102) 



and which satisfies the normalization condition, 



(103) 
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For simplicity, from here on we remove the initial-condition terms from our equations 
by assuming that cp^(O) = (p^(0), which implies that (p^ = (p^ holds at all times. 

We now move toward a Schmidt decomposition by separating out the transverse 
degrees of freedom from the other ones in Eq. (101), 



where 



^12 = Co(Po + eCi(pi, (104) 

Co = eoe-^^°^/^ (105) 
Ci = -e-^^°^/^f;^.^^^, (106) 



^0 = 4 I "J h (107) 




and 



^i=^fl'"!i^'"!iSiV (108) 

We can formulate the two-term decomposition (104) because, as in the static 
case [ ], we can segregate the n dependence of the n > 1 corrections in a transverse 
piece, leaving a longitudinal-hyperfine piece that has no n dependence. The resulting 
decomposition is not yet in Schmidt form, for the functions are not orthonormal, but 
we can easily transform it in Schmidt form. We begin by noting that 

(CoKo) = 1, (109) 
(Co|Ci)=0, (110) 



CO 



3\2 



(<Po|(Po) = 1, (112) 

(<po|(Pi) = ^ {giivl + ~g22vl + 2^12^7^) , (113) 

((Pi I (Pi) = ^ {glivla + ~gWLG + (2^11 + gi2)gi2VL^ + (2^22 + gi2)gi2vf^) ■ (114) 

We now write 

1^12 = (Co + e((po|(Pi)Ci)(Po + eCi((Pi - ((po|(Pi)(Po), (115) 
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Figure 23. Concurrence (125) of the two-mode condensate state (120) as a function 
of time and atom number for = 2. The entanglement remains remarkably small 
even for relatively large N ^ indicating that the condensate is well approximated by a 
single Schmidt term. In comparison to Figs. 24 and 25, the entanglement between the 
transverse dimensions and the longitudinal-hyperfine degrees of freedom decreases as 
the potential becomes more homogeneous. 



which inspires us to define the foUowing functions, 

Xo = Co + e((po|(Pi)Ci, (116) 

$0 = <Po, (118) 

= _^izMiSM^. (119) 

At the order we are working, Xo and Xi a-^e orthonormal, whereas 4)0 and cjji are 
orthonormal by construction. Thus Eqs. (116)-(119) form the perturbative Schmidt 
basis of the two-mode condensate state (104). Therefore, we write 

iji2 = VXoXo^o + e \/AiXi$i, (120) 
where the squared Schmidt coefficients are given by 

Ao = 1, (121) 

Ai = (CiKi)(((Pi|(Pi>-|((Po|(Pi>r). (122) 
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Figure 24. Concurrence (125) of the two-mode condensate state (120) as a function of 
time and atom number for q = A. As for a harmonic trap, the entanglement between the 
transverse dimensions and the longitudinal-hyperfine degrees of freedom is remarkably 
small even for relatively large A/", indicating that the condensate is well approximated 
by a single Schmidt term. 

As in [ ], we use Wootters's concurrence to quantify the bipartite entanglement 
of state (104). The concurrence of a pure state of two qubits, \^ ab)) varies smoothly 
from for product states to 1 for maximaUy entangled states. From its definition, 
C = K^Asl^y ® terms of the Pauh matrix Gy and the complex conjugate 

of \^ ab)) it is easy to show that the concurrence for the perturbative condensate single- 
particle state (120) is given by 

C = 2Va^ = 2(Ci|Ci) (((Pi|(Pi) - |((Po|(Pi)r) , (123) 

where we use the perturbative coefficients Aq and Ai as given by Eqs. (121) and (122). 
In the case of a cigar-shaped BEC, 

where we use the polylogarithm function Lis(z) = Yl^=i^^/'^^' Thus, it follows that 
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Figure 25. Concurrence (125) of the two-mode condensate state (120) as a function of 
time and atom number for q = In comparison to the other traps, the concurrence for 
(7 = 10 has the lowest values. Therefore, to very good approximation the nonseparable 
corrections are indeed negligible. 



the time-dependent concurrence can be written as 



, o~ f~ 112 , ~ 221 ~ 1 12 ~ 2 12\ 



(125) 



We show the concurrence (125) in Figs. 23, 24, and 25 as a function of time and 
atom number, for g = 2, 4, and 10. In all cases, the concurrence oscillates in time as well 
as increases as gets larger. More importantly, we can verify that the entanglement 
between the transverse and the remaining degrees of freedom is quite small for all cases, 
even for larger atom numbers. Moreover, it also decreases as the trap gets harder. 
These results indicate that the two-mode condensate state (120) is well described by 
the product state corresponding to the first Schmidt term. As we discussed before, the 
dominant term in the Schmidt decomposition corresponds to the optimal product-state 
approximation to state (100). This result not only confirms that the tightly confined 
dimensions can indeed be decoupled from the problem, but it naturally shows the 
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optimal way to do it. 
7. Conclusion 

We have presented a detailed theoretical analysis of the mean-field dynamics and the 
reduced-dimension character of two-mode BECs in highly anisotropic traps. These 
systems are promising candidates for the implementation of high precision nonlinear 
interferometers and our analysis shows the feasibility of a recent experimental proposal 
of a Ramsey interferometry scheme [ ] . 

In spite of the strong anisotropy of the condensate, our results also show that 
the eflFects of the nonlinear scattering interaction on the transverse degrees of freedom 
of the gas cannot be in general neglected. This brings into question the accuracy of 
the usual reduced-dimension (quasi- ID) approximation, both spatially and temporally. 
We have formulated a model that takes into account such 3D-induced effects by 
means of perturbative techniques. Because of the reshaping of the condensate in the 
transverse direction(s), additional eflFective three-body, attractive, intra- and inter-mode 
interactions have to be included in the reduced-dimension equations of the longitudinal 
dimension(s). Moreover, we have shown that these equations in fact provide the 
optimal reduced-dimension evolution for the two-mode condensate in general quasi- 
ID and quasi-2D geometries. We have demonstrated this result by developing a 
Schmidt decomposition of the two-mode condensate single-particle wave function, which 
separates out the transverse from the longitudinal-internal degrees of freedom. Using 
this formalism, we have verified that the entanglement between the transverse and 
the remaining degrees of freedom is quite small, even for larger atom numbers. In 
other words, this result shows that the two-mode condensate single-particle state is well 
described by the product state corresponding to the dominant term in the Schmidt 
decomposition. This corresponds to the optimal product-state approximation [ ]. 

Finally, as far as the interferometry process goes, our perturbative model takes into 
account the spatial diflFerentiation of the wave functions of the two modes and gives a 
significantly more accurate analytical account of the fringe signal than does the reduced- 
dimension approximation. We point out that this theoretical model is not limited to 
the study of interference eflFects, but can also be applied to the study of more general 
dynamics and diflFerent phenomena, such as soliton propagation. 

Acknowledgment s 

We thank A. Shaji for discussions at the beginning of this project. This work was 
supported in part by National Science Foundation Crant Nos. PHY-0903953 and PHY- 
1005540. 

References 

[1] A. M. Rey, L. Jiang, and M. D. Lukin, Phys. Rev. A 76, 053617 (2007). 



Mean-field dynamics of two-mode BECs in highly anisotropic potentials 



45 



[2] S. Boixo, A. Datta, M. J. Davis, S. T. Flammia, A. Shaji, and C. M. Caves, Phys. Rev. Lett. 

101, 040403 (2008). 
[3] S. Choi and B. Sundaram, Phys. Rev. A 77, 053613 (2008). 

[4] S. Boixo, A. Datta, M. J. Davis, A. Shaji, A. B. Tacla, and C. M. Caves, Phys. Rev. A 80, 032103 
(2009). 

[5] A. B. Tacla, S. Boixo, A. Datta, A. Shaji, and C. M. Caves, Phys. Rev. A 82, 053636 (2010). 

[6] A. B. Tacla, PhD Thesis, University of New Mexico (2012). 

[7] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006). 

[8] A. B. Tacla and C. M. Caves, Phys. Rev. A 84, 053606 (2011). 

[9] A. E. Muryshev, G. V. Shlyapnikov, W. Ertmer, K. Sengstock, and M. Lewenstein, Phys. Rev. 
Lett. 89, 110401 (2002). 
[10] S. Sinha, A. Y. Cherny, D. Kovrizhin, and J. Brand, Phys. Rev. Lett. 96, 030406 (2006). 
[11] L. Khaykovich and B. A. Malomed, Phys. Rev. A 74, 023607 (2006). 

[12] M. R. Matthews, D. S. Hah, D. S. Jin, J. R. Ensher, C. E. Wieman, E. A. Corneh, F. Dalfovo, 

C. Minniti, and S. Stringari, Phys. Rev. Lett. 81, 243 (1998). 

[13] D. S. Hah, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Corneh, Phys. Rev. Lett. 
81, 1539 (1998). 

[14] D. S. Hah, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Corneh, Phys. Rev. Lett. 
81, 1543 (1998). 

[15] K. M. Mertes, J. W. Merrill, R. Carretero-Gonzalez, D. J. Frantzeskakis, P. G. Kevrekidis, and 

D. S. Hah, Phys. Rev. Lett. 99, 190402 (2007). 

[16] R. P. Anderson, C. Ticknor, A. L Sidorov, and B. V. Hah, Phys. Rev. A 80, 023603 (2009). 
[17] M. Egorov, R. P. Anderson, V. Ivannikov, B. Opanchuk, P. Drummond, B. V. Hall, and A. L 

Sidorov, Phys. Rev. A 84, 021605(R) (2011). 
[18] M. Egorov, B. Opanchuk, P. Drummond, B. V. Hall, P. Hannaford, and A. L Sidorov, 

arXiv:1204.1591 (2012). 

[19] All numerical integrations are performed using the code generator XmdS^ available at 
www . xmds . org, except for the integration of the time-independent, three-dimensional CP 
equation. For this integration we use the method described in C. M. Dion and E. Cances, 
Comp. Phys. Comm. 177, 787 (2007). 

[20] In our simulations we assume that one can neglect the evolution of the spatial modes during the 
Ramsey pulses. According to [15, 16, 12], the two hyperfine levels of ^^Rb can be coupled by two- 
photon Rabi pulses whose duration is less than 1 ms, which is much faster than the characteristic 
time for the condensate to change its shape, thereby justifying such an approximation. 

[21] Plots that are in harmonic trap units have (i) transverse (longitudinal) lengths in units of po {zq 
for (7 = 2), (ii) transverse (longitudinal) marginal distributions in units of l/pg (l/^o foi" ^ = 2), 
and (iii) density in units of 1/{pqZo). 

[22] R. B. Lockhart and M. J. Steiner, Phys. Rev. A 65, 022107 (2002). 

[23] Nielsen M. and Chuang L, Quantum Computation and Quantum Information (Cambridge 

University Press, 2000). 
[24] W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998). 



